Intro to BioEnergeticFoodWebs
by Chris Griffiths, Eva Delmas and Andrew Beckerman, Dec. 2020.
This document follows on from 'Julia in VS Code #1, #2 and #3' and assumes that you're still working in your active project.
This document introduces the BioEnergeticFoodWebs.jl and EcologicalNetworks.jl packages. It demonstrates how to run the BioEnergetic Food Web (BEFW) model, how to vary variables of interest (e.g., productivity) and construct experiments designed to investigate the effect of different variables on population and community dynamics.
For those that are unfamilar with the BEFW and it's application in Julia, we advise checking out the MEE paper before we start. Remember, the BEFW model is also based on a system of differential equations and is solved using the same engine as the DifferentialEquations.jl package.
Packages
First, import your package manager, activate and instantiate (only neccessary if you've closed VS Code between tutorials):
x
begin import Pkg Pkg.activate(".") Pkg.instantiate()endThen install the BioEnergeticFoodWebs.jl, EcologicalNetworks.jl and JLD2.jl packages and let Julia know which packages you want to use:
x
begin # install Pkg.add("BioEnergeticFoodWebs") Pkg.add("EcologicalNetworks") Pkg.add("JLD2") Pkg.add("Statistics") # use using BioEnergeticFoodWebs using EcologicalNetworks using JLD2 using Statistics using Plots using CSV using DataFrames using RandomendThe JLD2.jl package will be useful later as it allows you to directly export and load a BEFW output object. Let's also set a random seed for reproducibility:
0x00000015
163718196
1073255743
-1451494363
1072899844
-1618160924
1073258805
-1517321307
1073339839
-1150166203
1072942963
-1583246675
1072950051
-407139083
1073171397
1426608413
1073390650
507415766
1073502426
-1219197164
1073024304
-1705970718
1072820091
394465429
1073380663
887692544
1617232906
-1051299616
1513956225
382
0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
0x00000000000000000000000000000000
1002
0
xxxxxxxxxxRandom.seed!(21)Preamble
One of main advantages of running food web models in Julia is that simulations are fast and can be readily stored in your active project. With this in mind, make a new folder in your project called out_objects (right click > New Folder). Alternatively, you can create an out_objects folder directly using mkdir().
Running the BEFW
There are four major steps when running the BioEnergetic Food Web model in Julia:
Generate an initial network
Fix parameters
Simulate
Explore output and plot
Initial network
Before running the BEFW model, we have to construct an initial random network using the niche model. The network is characterised by the number of species in the network and its connectance value. Here, we generate a network of 20 species with a connectance value of 0.15:
x
md"## Running the BEFWThere are four major steps when running the BioEnergetic Food Web model in Julia:1. Generate an initial network 2. Fix parameters3. Simulate4. Explore output and plot### Initial network Before running the BEFW model, we have to construct an initial random network using [the niche model](https://www.nature.com/articles/35004572?cacheBust=1510239451067). The network is characterised by the number of species in the network and its [connectance](https://en.wikipedia.org/wiki/Ecological_network) value. Here, we generate a network of 20 species with a connectance value of 0.15:"20×20 Array{Int64,2}:
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
⋮ ⋮ ⋮ ⋮
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0
0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1
0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0x
begin # generate network A_bool = EcologicalNetworks.nichemodel(20,0.15) # convert the UnipartiteNetwork object into a matrix of 1s and 0s A = Int.(A_bool.A) # 1s indicate an interaction among species and 0s no interaction. In the packages used here, the networks are directed from i to j (i eats j), describing the direction of the interaction, not of the flow of biomass.endYou can check the connectance of A using:
0.145x
# calculate connectanceco = sum(A)/(size(A,1)^2)Parameters
Prior to running the BEFW model, you have to create a vector of model parameters using the model_parameters function. Numerous parameter values can be specified within the model_parameters function, however, most of them have default values that are built into the BioEnergeticFoodWebs.jl package. For simplicity, we use the default values here:
:α
1.0
:e_carnivore
0.85
:Γh
0.0
0.0
0.0
0.5
0.0
0.5
0.5
0.5
0.5
0.0
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
:m_producer
1.0
:c
0.0
:h
1.0
:vertebrates
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
:tmpA
:rewire_method
:none
:trophic_rank
1.0
1.0
1.0
2.0
1.0
2.2
2.75
2.5
2.5
1.0
3.0
3.0
4.0
2.5
3.5
3.0
2.75
2.5
3.5
2.5
:w
20×20 Array{Float64,2}:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.2 0.2 0.2 0.2 0.2 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.25 0.25 0.25 0.0 0.0 0.0 0.0 0.0
⋮ ⋮ ⋱ ⋮
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.125 0.125 0.125 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.166667 … 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.25 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.125 0.125 0.125 0.125 0.125
0.0 0.0 0.25 0.25 0.25 0.25 0.0 0.0 0.0 0.0 0.0 0.0:TSR_type
:no_response
:e_herbivore
0.45
:productivity
:species
:efficiency
20×20 Array{Float64,2}:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.45 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.45 0.45 0.45 0.85 0.45 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.45 0.85 0.85 0.0 0.0 0.0 0.0 0.0 0.0
⋮ ⋮ ⋱ ⋮
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.85 0.85 0.85 0.85 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.85 … 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.85 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.85 0.85 0.85 0.85 0.85 0.85
0.0 0.0 0.45 0.85 0.45 0.85 0.0 0.0 0.0 0.0 0.0 0.0 0.0:K
1.0
:S
20
:extinctionstime
:is_producer
true
true
true
false
true
false
false
false
false
true
false
false
false
false
false
false
false
false
false
false
:dp
#10
:x
0.138
0.138
0.138
0.3141
0.138
0.3141
0.3141
0.3141
0.3141
0.138
0.3141
0.3141
0.3141
0.3141
0.3141
0.3141
0.3141
0.3141
0.3141
0.3141
:Z
1.0
:extinctions
:dry_mass_293
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
:np
5
:A
20×20 Array{Int64,2}:
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
⋮ ⋮ ⋮ ⋮
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0
0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1
0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0:Γ
0.0
0.0
0.0
0.5
0.0
0.5
0.5
0.5
0.5
0.0
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
:ar
0.0
0.0
0.0
16.0
0.0
16.0
16.0
16.0
16.0
0.0
16.0
16.0
16.0
16.0
16.0
16.0
16.0
16.0
16.0
16.0
:bodymass
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
:is_herbivore
false
false
false
true
false
true
true
true
true
false
false
true
false
true
false
true
true
true
false
true
x
# create model parametersp = model_parameters(A)# in the most simple case, the model_parameters function simply requires AFor more information and a full list of the parameters and their defaults values type ?model_parameters in the REPL.
Simulate
To run the BEFW model, we first assign biomasses at random to each species and then simulate the biomass dynamics forward using the simulate function:
:p
:α
1.0
:e_carnivore
0.85
:Γh
0.0
0.0
0.0
0.5
0.0
0.5
0.5
0.5
0.5
0.0
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
:m_producer
1.0
:c
0.0
:h
1.0
:vertebrates
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
false
:tmpA
20×20 Array{Int64,2}:
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
⋮ ⋮ ⋮ ⋮
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0
0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1
0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 020×20 Array{Int64,2}:
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
⋮ ⋮ ⋮ ⋮
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0
0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1
0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 020×20 Array{Int64,2}:
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
⋮ ⋮ ⋮ ⋮
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0
0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1
0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 020×20 Array{Int64,2}:
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
⋮ ⋮ ⋮ ⋮
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0
0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1
0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 020×20 Array{Int64,2}:
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
⋮ ⋮ ⋮ ⋮
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0
0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1
0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 020×20 Array{Int64,2}:
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
⋮ ⋮ ⋮ ⋮
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0
0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1
0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 020×20 Array{Int64,2}:
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
⋮ ⋮ ⋮ ⋮
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0
0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1
0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 020×20 Array{Int64,2}:
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
⋮ ⋮ ⋮ ⋮
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0
0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1
0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 020×20 Array{Int64,2}:
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
⋮ ⋮ ⋮ ⋮
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0
0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1
0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0:rewire_method
:none
:trophic_rank
1.0
1.0
1.0
2.0
1.0
2.2
2.75
2.5
2.5
1.0
3.0
3.0
4.0
2.5
3.5
3.0
2.75
2.5
3.5
2.5
:w
20×20 Array{Float64,2}:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.2 0.2 0.2 0.2 0.2 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.25 0.25 0.25 0.0 0.0 0.0 0.0 0.0
⋮ ⋮ ⋱ ⋮
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.125 0.125 0.125 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.166667 … 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.25 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.125 0.125 0.125 0.125 0.125
0.0 0.0 0.25 0.25 0.25 0.25 0.0 0.0 0.0 0.0 0.0 0.0:TSR_type
:no_response
:e_herbivore
0.45
:productivity
:species
:efficiency
20×20 Array{Float64,2}:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.45 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.45 0.45 0.45 0.85 0.45 0.0 0.0 … 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.45 0.85 0.85 0.0 0.0 0.0 0.0 0.0 0.0
⋮ ⋮ ⋱ ⋮
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.85 0.85 0.85 0.85 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.85 … 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.85 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.85 0.85 0.85 0.85 0.85 0.85
0.0 0.0 0.45 0.85 0.45 0.85 0.0 0.0 0.0 0.0 0.0 0.0 0.0:K
1.0
:S
20
:extinctionstime
11.0425
11
19.7522
4
29.6163
13
45.5946
15
76.8489
16
98.1864
5
98.1864
12
134.606
17
205.609
19
:is_producer
true
true
true
false
true
false
false
false
false
true
false
false
false
false
false
false
false
false
false
false
:dp
#10
:x
0.138
0.138
0.138
0.3141
0.138
0.3141
0.3141
0.3141
0.3141
0.138
0.3141
0.3141
0.3141
0.3141
0.3141
0.3141
0.3141
0.3141
0.3141
0.3141
:Z
1.0
:extinctions
4
5
11
12
13
15
16
17
19
:dry_mass_293
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
0.153846
:np
5
:A
20×20 Array{Int64,2}:
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
⋮ ⋮ ⋮ ⋮
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0
0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0
0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1
0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0:Γ
0.0
0.0
0.0
0.5
0.0
0.5
0.5
0.5
0.5
0.0
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
:ar
0.0
0.0
0.0
16.0
0.0
16.0
16.0
16.0
16.0
0.0
16.0
16.0
16.0
16.0
16.0
16.0
16.0
16.0
16.0
16.0
:bodymass
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
:is_herbivore
false
false
false
true
false
true
true
true
true
false
false
true
false
true
false
true
true
true
false
true
:B
8017×20 Adjoint{Float64,Array{Float64,2}}:
0.524056 0.535102 0.977416 0.906369 … 0.459164 0.615476 0.245534 0.490372
0.517543 0.527093 0.192875 0.923082 0.553908 0.718219 0.293288 0.625466
0.50712 0.515335 0.0200148 0.710628 0.614887 0.701168 0.358189 0.746566
0.501974 0.509129 0.00282734 0.491415 0.639262 0.621707 0.440936 0.839793
0.505704 0.51205 0.000635531 0.328812 0.626891 0.530129 0.540591 0.893063
0.518701 0.524421 0.000204033 0.217314 … 0.581946 0.441034 0.654733 0.901809
0.539661 0.544868 8.4954e-5 0.144 0.51424 0.358944 0.77923 0.870207
⋮ ⋱
0.559053 0.559053 1.60899e-5 0.0 0.0 0.139536 0.0 0.356004
0.559054 0.559054 1.54769e-5 0.0 0.0 0.139536 0.0 0.356005
0.559054 0.559054 1.45103e-5 0.0 0.0 0.139536 0.0 0.356006
0.559055 0.559055 1.33312e-5 0.0 0.0 0.139536 0.0 0.356007
0.559055 0.559055 1.21177e-5 0.0 … 0.0 0.139536 0.0 0.356008
0.559055 0.559055 1.10852e-5 0.0 0.0 0.139536 0.0 0.356008:t
0.0
0.25
0.5
0.75
1.0
1.25
1.5
1.75
2.0
2.25
2.5
2.75
3.0
3.25
3.5
3.75
4.0
4.25
4.5
4.75
1997.75
1998.0
1998.25
1998.5
1998.75
1999.0
1999.25
1999.5
1999.75
2000.0
xxxxxxxxxxbegin # assign biomasses bm = rand(size(A,1)) # select biomasses at random between ]0:1[ # simulate out = simulate(p, bm, start=0, stop=2000) # this might take a few secondsendThe simulate function requires the model parameters p and the species biomasses bm. In addition, you can specify the timespan of the simulation (using the start and stop arguments), fix a species extinction threshold (using extinction_threshold) and select a solver (using use). For more information type ?simulate in the REPL.
Output and plot
Once the simulation finishes, the output is stored as a dictionary called out. Within out there are three entries:
out[:p]- lists the parametersout[:B]- biomass of each species through timeout[:t]- timesteps (these typically increase in 0.25 intervals)
The biomass dynamics of each species can then be plotted. Similar to the DifferentialEquations.jl package, the BioEnergeticFoodWebs.jl package also has it's own built in plotting recipe:
x
# plotPlots.plot(out[:t], out[:B], legend = true, ylabel = "Biomass", xlabel = "Time")# this may take a minute to renderYou'll notice that the biomass dynamics are noisey during the first few hundred time steps, these are the system's transient dynamics. The dynamics then settle into a steady state where the system can be assumed to be at equilbirum. You'll also notice that some species go extinct and some persist, the number of species in the food web can found using out[:p][:S] and the identity of those that went extinct using out[:p][:extinctions].
The BioEnergeticFoodWebs.jl package also has a range of built in functions that conveniently calculate some of the key metrics of the food web, these include the total biomass, the diversity, the species persistence and the temporal stability:
2.0796559291208965x
# total biomassbiomass = total_biomass(out, last=1000)0.7182050496471224xxxxxxxxxx# diversitydiversity = foodweb_evenness(out, last=1000)0.55xxxxxxxxxx# persistencepersistence = species_persistence(out, last=1000)-0.04410056440077254xxxxxxxxxx# stability stability = population_stability(out, last=1000)Each of these functions will output a single value. This value is the average over the last 1000 time steps. For more information, use ? to access the help files on each function in the REPL (e.g., ?species_persistence).
Variables
Once you've got the BEFW model running, the next step is to vary a variable of interest and rerun. For example, we might be interested in what affect a small change in Z (consumer-resource body mass ratio) has on the estimated food web and its biomass dynamics. The default value for Z is 1.0, but what happens if we increase it to 10.0:
xxxxxxxxxxbegin # set Z Z = 10.0 # create model parameters p_z = model_parameters(A, Z = Z) # assign biomasses bm_z = rand(size(A,1)) # simulate out_z = simulate(p_z, bm_z, start=0, stop=2000) # plot Plots.plot(out_z[:t], out_z[:B], legend = true, ylabel = "Biomass", xlabel = "Time")endSimilarly, what happens if we also increase the carrying capacity (K) of the resource from 1.0 (default) to 5.0:
xxxxxxxxxxbegin # set K K = 5.0 # create model parameters p_K = model_parameters(A, Z = Z, K = K) # assign biomasses bm_K = rand(size(A,1)) # simulate out_K = simulate(p_K, bm_K, start=0, stop=2000) # plot Plots.plot(out_K[:t], out_K[:B], legend = true, ylabel = "Biomass", xlabel = "Time")endAs you've probably guessed, the main message here is that many variables can be changed in the BEFW model and it's super easy to do so. In the next step, we take this one step further.
Experiments
The next step is to construct a computional experiment designed to investigate the effect of different variables on population and community dynamics. To do this we construct a gradient of variables as vectors and then simulate the BEFW model multiple times using a loop. To illustrate this, we're going to reproduce example 1 from Delmas et al. 2016. The aim of this example is to investigate the effect of increasing K on food web diversity. In addition, we're also going to allow α (interspecific competition relative to intraspecific competition) to vary and repeat the experiment 5 times with 5 different initial networks.
First, we define the experiment by creating vectors of our variables and fixing the number of repetitions:
0.92
1.0
1.08
x
α = [0.92, 1.0, 1.08]# 0.92 - the interspecific competition is smaller than the intraspecific competition promoting coexistence# 1.0 - neutrally stable # 1.08 - the intraspecific competition is smaller the interspecific competition favouring competitive exclusion0.1
0.16681
0.278256
0.464159
0.774264
1.29155
2.15443
3.59381
5.99484
10.0
x
# vector of Kk = exp10.(range(-1,1,length=10))# log scale from 0.1 to 105x
# number of repsreps = 5We then create a dataframe to store the outputs:
| α | K | network | diversity | stability | biomass | |
|---|---|---|---|---|---|---|
| Any | Any | Any | Any | Any | Any |
xxxxxxxxxx# dataframedf = DataFrame(α = [], K = [], network = [], diversity = [], stability = [], biomass = [])and construct a while loop to generate the 5 unique initial networks, each of which contains 20 species with a connectance value of 0.15:
xxxxxxxxxxbegin # list to store networks global networks = [] # monitoring variable global l = length(networks) # while loop while l < reps # generate network A_bool = EcologicalNetworks.nichemodel(20,0.15) # convert the UnipartiteNetwork object into a matrix of 1s and 0s A = Int.(A_bool.A) # calculate connectance co = sum(A)/(size(A,1)^2) # ensure that connectance = 0.15 if co == 0.15 push!(networks, A) # save network is co = 0.15 end global l = length(networks) endendWe can then run the simulations by looping, using nested for loops, over the unique values of α and K, as well as the 5 unique initial networks. After each simulation we will save each output object to our active project as a JLD2 file and store any output metrics of interest in our dataframe:
xxxxxxxxxx# loop over networksfor h in 1:reps A = networks[h] # here, you might want to save a copy of the initial network using writedlm(A) # loop over α for i in 1:length(α) # loop over K for j in 1:length(k) # create model parameters p = model_parameters(A, α = α[i], K = k[j]) # assign biomasses bm = rand(size(A,1)) # simulate out = simulate(p, bm, start=0, stop=2000) # dummy naming variables α_num = α[i] K_num = k[j] # save `out` as a JLD2 object using the @save macro: "out_objects/model_output, network = $h, alpha = $α_num, K = $K_num.jld2" out # calculate output metrics diversity = foodweb_evenness(out, last = 1000) stability = population_stability(out, last = 1000) biomass = total_biomass(out, last = 1000) # push to df push!(df, [α[i], k[j], h, diversity, stability, biomass]) # print some stuff - see how the simulation is progressing println(("α = $α_num", "K = $K_num", "network = $h")) end endend# the code will be much faster if you remove the @save commandWe can then explore the outputs and plot our results. Here, instead of using the built in plotting recipe, we will construct a plot that matches figure 1 in Delmas et al. 2016. Specifically, we will plot food web diversity (y-axis) as a function of K (x-axis) and α (colour):
| variable | mean | min | median | max | nunique | nmissing | eltype | |
|---|---|---|---|---|---|---|---|---|
| Symbol | Float64 | Float64 | Float64 | Int64 | Int64 | DataType | ||
| 1 | :α | 1.0 | 0.92 | 1.0 | 1.08 | 3 | 0 | Any |
| 2 | :K | 5.0 | 5.0 | 5.0 | 5.0 | 1 | 0 | Any |
| 3 | :network | 3.0 | 1.0 | 3.0 | 5.0 | 5 | 0 | Any |
| 4 | :diversity | NaN | NaN | nothing | NaN | 9 | 0 | Any |
| 5 | :stability | -0.099386 | -1.3342 | 0.0 | -0.0 | 8 | 0 | Any |
| 6 | :biomass | 6.96612 | 2.99178 | 5.0 | 14.2182 | 9 | 0 | Any |
xxxxxxxxxx# explore outputdescribe(df)| α | K | network | diversity | stability | biomass | |
|---|---|---|---|---|---|---|
| Any | Any | Any | Any | Any | Any | |
| 1 | 0.92 | 5.0 | 1.0 | NaN | -0.0 | 5.0 |
| 2 | 1.0 | 5.0 | 1.0 | NaN | -0.0 | 5.0 |
| 3 | 1.08 | 5.0 | 1.0 | NaN | -0.0 | 5.0 |
| 4 | 0.92 | 5.0 | 2.0 | 0.928022 | -0.00466869 | 2.99184 |
| 5 | 1.0 | 5.0 | 2.0 | 0.948159 | -0.00468039 | 2.99178 |
| 6 | 1.08 | 5.0 | 2.0 | 0.977442 | -0.00468403 | 2.99182 |
xxxxxxxxxxfirst(df,6)| α | K | network | diversity | stability | biomass | |
|---|---|---|---|---|---|---|
| Any | Any | Any | Any | Any | Any | |
| 1 | 0.92 | 5.0 | 4.0 | NaN | -0.0 | 5.0 |
| 2 | 1.0 | 5.0 | 4.0 | NaN | -0.0 | 5.0 |
| 3 | 1.08 | 5.0 | 4.0 | 0.413271 | -1.3342 | 7.86188 |
| 4 | 0.92 | 5.0 | 5.0 | NaN | -0.0 | 5.0 |
| 5 | 1.0 | 5.0 | 5.0 | NaN | -0.0 | 5.0 |
| 6 | 1.08 | 5.0 | 5.0 | 1.0 | -0.0 | 10.0 |
xxxxxxxxxxlast(df,6)xxxxxxxxxx# plot# initialise an empty plotpl = Plots.plot([NaN], [NaN], label = "", ylims = (0,1.1), leg = :bottomright, foreground_colour_legend = nothing, xticks = (log10.(k), string.(round.(k, digits = 1))), xlabel = "Carrying capacity", ylabel = "Food web diversity (evenness)"):square
:diamond
:utriangle
xxxxxxxxxx# set marker shapesshp = [:square, :diamond, :utriangle]:solid
:dash
:dot
xxxxxxxxxx# set line typesls = [:solid, :dash, :dot]
xxxxxxxxxx# set coloursclr = [RGB(174/255, 139/255, 194/255), RGB(188/255, 188/255, 188/255), RGB(124/255, 189/255, 122/255)]# when we define colours in Julia they are printed "Coexistence"
"Neutral"
"Exclusion"
x
# set legend labelslbl = ["Coexistence", "Neutral", "Exclusion"]xxxxxxxxxx# make the plotfor (i, α) in enumerate(α) # subset tmp = df[df.α .== α, :] # remove NaN values tmp = tmp[.!(isnan.(tmp.diversity)), :] # calculate mean across reps meandf = by(tmp, :K, :diversity => mean) # command to avoid printing legends multiple times l = i == 1 ? lbl[i] : "" # add to pl plot!(pl, log10.(meandf.K), meandf.diversity_mean, msc = clr[i], mc = :white, msw = 3, markershape = shp[i], linestyle = ls[i], lc = clr[i], lw = 2, label = lbl[i], seriestype = [:line :scatter])endxxxxxxxxxx# display plotplot(pl)Finally, we can save our dataframe as a .csv file:
"My_data.csv"xxxxxxxxxx# saveCSV.write("My_data.csv", df)